USING  LINEAR  PROGRAMMING  AS  A  SIMPLEX  SUBROUTINE 

f 


R.  J.  Clasen 

November  1965 


h  ■ 

P-3267 


Approved  for  release  Py  the  Clearinghouse  lor 
Federal  Scientific  and  Technical  Information 


ABSTRACT 


This  Paper  discusses  the  problems  involved  in  using 
linear  programming  as  a  subroutine  of  a  larger  routine. 
Proposals  are  made  for  eliminating  the  tolerance  selection 
problem,  and  for  improving  the  accuracy  of  inversions. 

Sample  programs  are  given  in  FORTRAN  IV  and  in  ALGOL. 


-1- 


USING  LINEAR  PROGRAMMING  AS  A  SIMPLEX  SUBROUTINE 

R.  J.  Clasen 


The  RAND  Corporation,  Santa  Monica,  California 


INTRODUCTION 


Use  of  linear  programming  has  increased  rapidly  in 
recent  years.  The  major  use  was  and  still  is  to  solve  a 
problem  stated  in  terms  of  linear  programming;  generally, 
only  the  answer  to  the  problem  was  necessary.  The  one 
difficulty,  usually,  was  preparing  the  input  in  a  suitable 
format,  and  obtaining  the  output  answer  in  an  equally 
suitable  format. 

However,  a  problem  occasionally  arises  of  which  only 
a  small  part  involves  obtaining  the  answer  to  a  linear 
programming  problem.  For  example,  one  may  use  linear  pro¬ 
gramming  to  initially  estimate  the  solution  to  a  chemical 
equilibrium  problem  rl\  In  this  case,  it  would  be  con¬ 
venient  to  have  linear  programming  as  a  subroutine.  We 


+ 
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define  our  linear  programming  problem  as  that  of  deter 

mining  x, such  that 
°  12  n 


c  .x. 
J  J 


j  =  l 


is  a  minimum  over  all  sets  x, ,x„,...,x  that  satisfy 

12  n 


n 

Y,  aijXj  =  bi  •  •  •  »m  » 

j-1 

and 


Xj  ^  0  j-1 , 2 , . . . , n  . 

The  details  of  the  construction  of  such  a  subroutine  are 
the  subject  of  this  Paper.  The  criteria  for  determining 
what  procedures  should  be  used  are: 

1)  Accuracy  of  the  solution,  including  its  inde¬ 
pendence  of  row  and  column  scale  factors, 

2)  Storage  needed  for  data, 

3)  Length  of  the  subroutine  needed, 

4)  Speed  of  solution-- 

in  that  order. 

DETERMINING  TOLERANCES 

A  person  solving  linear  programming  problems  with  the 
simplex  method  "2  generally  has  trouble  estimating  certain 
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tolerances.  Typically,  the  following  three  tolerances  are 
needed : 

1)  The  Pivot  Tolerance:  a  number  in  the  pivot  column 
is  considered  zero  if  its  absolute  value  does  not 
exceed  the  pivot  tolerance. 

2)  The  Zero  Tolerance:  a  number  in  the  solut ion 
vector  is  considered  zero  if  its  absolute  value 
does  not  exceed  the  zero  tolerance. 

3)  The  Cost  Tolerance:  a  reduced  cost  is  considered 
zero  if  its  absolute  value  does  not  exceed  the 
cost  tolerance. 


We  abbreviate  these  tolerances  as  TP,  TZ ,  and  TC , 
respectively . 

When  using  a  simplex  subroutine,  the  user  is  unable 
to  see  his  input  numbers  before  he  uses  the  subroutine; 
hence,  he  is  not  able  to  determine  the  proper  tolerances. 
Therefore,  the  simplex  subroutine  should  either  calculate 
its  own  tolerances,  or  use  a  method  that  does  not  need 
tolerances.  We  used  a  combination  of  these  two  procedures. 

Before  giving  the  procedure  for  obtaining  tolerances, 
let  us  first  review  the  basic  steps  in  the  simplex  method 
using  the  "explicit  inverse"  method: 

1)  Check  tne  solution  vector  for  feasibility.  Let 
the  basic  part  of  the  solution  vector  be  denoted 


This  is  what  Cutler  and  Wolfe  [3]  call  x.  ,x.  , 


••• »xj 


h  J2 


m 
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2)  Calculate  the  prices--the  "phase  one"  prices  if 
the  problem  is  not  yet  feasible,  the  "phase  two" 
prices  if  the  problem  is  feasible. 

3)  Calculate  the  reduced  costs  and  find  the  column, 
JT,  with  the  minimum  reduced  cost,  MRC.  If 

MRC  0,  terminate  the  subroutine;  if  MRC  0, 
the  pivot  column  is  column  JT. 

4)  Obtain  the  column  vector  JT  by  multiplying  the 

inverse  and  the  original  column  JT.  Let  the 

column  obtained  be  called  y,,y0>---*y  • 

12  m 

3)  Obtain  the  pivot  row,  IR,  by  using  the  subscript 

w. 

that  causes  the  quantity  —  to  be  a  minimum  for 

w . 

all  non-zero  y.  for  which  —  0.  Slight  varia- 

L  y  i 

tions  of  this  rule  may  be  used  when  w^  ^  0.  If 

no  row  is  found,  then  for  each  i  either  y.  =0 

i 

w . 

or  —  ^  0.  In  this  case,  we  have  an  infinite 

solution  and  the  subroutine  is  terminated. 

6)  Update  the  inverse,  the  "phase  two"  prices,  and 
the  w^  by  pivoting  on  (IR,  JT). 

These  steps  are  typically  repeated  until  the  subroutine 

terminates  at  either  Step  3  or  Step  5.  The  initial  basis 

may  be  vacuous  (or  "artificial")  and  the  initial  inverse 

may  be  the  identity.  In  addition  to  these  steps,  every 

m/ 2  to  m  iterations,  it  is  usually  desirable  to  "re-invert" 

the  basis,  so  that  the  round-off  error  is  not  too  large. 

If  the  re- inversion  is  to  be  done  every  NVER  times,  we 

adjoin  to  this  system  a  counter,  INVC  (which  is  zero 

initially),  and  the  following  step: 
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7)  Increase  IN7C  by  1.  If  INVC  <  NVER ,  go  to  Step  1. 
Otherwise  set  INVC  to  zero  and  invert  the  basis, 
then  go  to  Step  1. 

The  pivot  tolerance  (TP)  is  used  in  Step  5  to  determine 

whether  or  not  a  is  zero.  A  satisfactory  method  of 

computing  this  tolerance  is  first  to  compute,  on  every 

m 

iteration,  YMAX  =  max  y. I.  Then  the  pivot  tolerance  for 

i=l  1 

that  iteration  is  taken  to  be  TP  =  YMAX*2  .  Then,  on  a 

machine  that  carries  numbers  in  floating  point  notation  to 

-27 

a  relative  accuracy  of  2  ,  we  assume  y^  =  0  in  Step  5  if 

| y ^ !  £  TP.  Thus  we  are,  in  effect,  assuming  the  last  11 
bits  are  not  significant.  While  this  is  an  ad  hoc  rule, 
it  has  worked  reasonably  well  for  problems  on  the  order 
of  50  constraints. 

The  zero  tolerance  (TZ)  is  typically  used  in  Step  1 
to  eliminate  small  negative  numbers  when  determining 
feasibility.  We  do  not  calculate  a  zero  tolerance,  but 
we  eliminate  the  small  negative  numbers  at  the  source, 
i.e.,  in  Step  6--the  pivot.  In  other  words,  when  we 
calculate  a  new  solution  vector  w|,  we  keep  the  old  solu¬ 
tion  vector  w then,  if  w.  ?  0  and  w!  0,  we  set  w!  to 

l  l  l  -  i 

zero.  These  negative  w^  can  also  arise  after  a  reinversion. 
If  the  reinversion  produces  a  full  basis,  and  the  problem 
was  feasible  before  the  reinversion  began,  we  set  to  zero 
any  negative  w^  generated  by  the  inversion. 
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Step  3  requires  the  cost  tolerance  (TC)  to  determine 
whether  a  small  negative  cost  is  "really"  zero.  "Really" 
zero  means  the  number  would  be  zero  if  infinite  computing 
accuracy  were  used.  The  present  method  is  to  temporarily 
ignore  the  possibility  that  a  reduced  cost  is  a  rounded 
zero,  and  to  find  the  minimum  reduced  cost  using  no  (or  a 
zero)  tolerance.  Then  we  proceed  as  usual  to  Step  4  to 
obtain  the  prospective  pivot  column  and  the  pivot  tolerance. 
We  use  the  number  TP  for  TC,  but  we  do  not  terminate  the 
first  time  the  minimum  reduced  cost  exceeds  -TC.  Instead, 
we  require  the  minimum  reduced  cost  to  exceed  -TC  for  two 
consecutive  iterations. 

Now,  if  the  column  to  be  pivoted  on  "really"  had  a 
zero  reduced  cost,  we  may  find  no  pivot  row  in  Step  3.  If 
we  do  find  a  pivot  row  in  a  column  with  a  reduced  cost  of 
zero,  it  does  no  harm  (except  possibly  some  lost  time)  to 
pivot  in  this  column,  since  we  will  merely  obtain  an  alter¬ 
nate  optimal  solution.  If  in  Step  5  we  find  no  pivot  row, 
we  do  not  declare  an  infinite  solution  unless  the  problem 
is  feasible,  and  the  minimum  reduced  cosr  (MRC)  calculated 
in  Step  3  satisfies:  -  1000*  MRC  <  YMAX.  If  these  con¬ 
ditions  are  not  satisfied,  we  declare  an  optimal  solution 
if  the  problem  is  in  a  feasible  state,  or  a  "no  feasible 
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solution"  if  the  problem  is  in  an  infeasible  state.  This 
additional  test  effectively  prevents  a  false  declaration 
of  an  "infinite  solution." 

Some  linear  programming  routines  allow  redundant  con¬ 
straints.  Because  we  use  an  effective  zero  tolerance  of 
zero,  we  do  not  allow  redundant  rows.  Hence  if  we  are 
not  able  to  pivot  in  every  row,  we  declare  an  infeasible 
solution.  Thus,  a  problem  is  infeasible  if  there  exist 
numbers  P^»P2,,,,’Pm  being  t^ie  number  of  constraints) 
not  all  zero  such  that 

m 

l  P^  *  0  for  j  =1 , 2 , 3 , .  .  .  ,n 
i=l 


and 


m 

l  pibi  2  0  • 

i-1 

If  we  have  no  feasible  solution  to  a  problem,  the  shadow 
prices  obtained  will  be  proportional  to  the  above  p^. 
(The  proportion  may  be  the  negative  of  the  p^  above,  in 
which  case  both  inequalities  above  would  be  reversed.) 
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INVERSION  METHOD 

The  method  of  inverting  in  our  subroutine  differs 
from  the  usual  procedure.  A  typical  method  of  inverting 
is  to  pivot  in  each  column  that  is  part  of  the  basis, 
pivoting  in  the  row  in  which  the  column  entry  has  the 
largest  absolute  value.  For  example,  suppose 


where  c  is  a  small  number  whose  absolute  value  is  less 

than  half  the  accuracy  to  which  unity  can  be  represented. 

-Q 

Thus,  on  an  IBM  7090,  e  might  be  10  .  Now  pivoting  in 

the  first  column,  we  would  chose  the  first  row  and  obtain: 


I 

1 

0 

f  \\ 

II 

l< 

,  B  = 

I 

^  0 

J 

But  we  assume  that  c  is  small  enough  so  that  e-^  numerically 
becomes  - \ .  Then  we  pivot  in  the  second  row  of  the  second 


column  and  obtain: 
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f  0  N 
l 1 ) 

hence,  =  0,  =  1.  We  note  the  correct  answer  is 

=  c,  w2  =  l-2(.  The  (  has  been  absorbed  into  the 
round-off  error. 

The  following  pivot  scheme  avoids  this  problem: 

1)  Let  yi,y2,. . . ,y  be  the  transformed  column 

in  which  we  wish  to  pivot.  Let  b,  ,b_, .  .  .  ,b  be 

r  12  m 

the  current  values  of  the  transformed  constant 

vector  (i.e.,  what  was  previously  called  w^) . 

2)  Let  S  be  the  set  of  integers  such  that  icS  if 
|y^  |  >  TP  where  TP  is  the  pivot  tolerance  calcu¬ 
lated  as  above. 

3)  Let  T  be  the  set  of  integers  such  that  icT  if 
b.  =  0. 

l 

A)  If  SOT  is  not  vacuous,  do  Step  A;  if  SPT  is 
vacuous,  do  Step  B. 


A)  Let  IR ,  IRc  SPT,  be  an  integer  such  that 


IR 


^  |y^l  for  all  ie  SDT. 


B)  Let  IR ,  IReS,  be  an  integer  such  that 

for  all  icS. 


V1R 

*1 

£ir 

6i 

—  sS** 
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Thus  we  are  effectively  pivoting  in  the  row  that  has  the 
largest  ratio  ’y^/b^ 1 .  In  the  above  example,  the  largest 
ratio  in  the  first  column  is  the  second  row;  hence  the  re¬ 
sult  of  the  first  pivot  is 


Here  the  l-2c  will  be  rounded,  presumably,  to  unity.  Then 
we  pivot  in  the  first  row  of  the  second  column  to  obtain: 


Hence  we  get  w^  =  1,  w^  =  f,  which  is  actually  the  correct 
answer  within  round-off  error.  This  method  has  the  dis¬ 
advantage  that,  in  not  choosing  the  largest  element  in  a 
column  to  determine  the  pivot  row,  we  may  get  an  element 
that  barely  exceeds  the  pivot  tolerance.  But  we  nonethe¬ 
less  believe  that  the  increased  accuracy  justifies  this 
approach. 
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Other  choices  could  have  been  made  for  the  pivot 

element.  For  example,  instead  of  pivoting  on  columns  in 

column  order,  one  could  first  choose  the  eligible  row 

with  the  smallest  w.,  then  choose  from  that,  row  of  the 

J 

matrix  the  elegible  element  with  the  largest  absolute 
value.  This  was  not  done  because  of  timing  considerations, 
and  because  it  would  make  the  selections  scale-dependent. 
Note  that,  in  our  method,  the  only  scale-dependent  opera¬ 
tion  is  finding  the  pivot  column  in  Step  3. 

4- 

The  simplex  subroutine  described  here  is  available 
as  a  FORTRAN  IV  routine  through  the  IBM-users  group  (SHARK). 
The  cards  for  this  subroutine  may  be  obtained  by  writing 
directly  to 

SHARK  Distribution  Agency 
DP  Program  Information  Department 
IBM  Corporation 
40  Saw  Mill  River  Road 
Hawthorne,  New  York  10532 

and  asking  for  SHARE  Distribution  Agency  No.  3384. 

Employees  of  RAND  may  obtain  a  copy  from  the  SHARI 
program  librarian.  Pearl  Leonhardt ,  by  asking  for  the 
W026  routine. 


Appendix 

THE  FORTRAN  SUBROUTINE 


The  simplex  subroutine,  SIMPLE,  may  be  used  to  solve 

a  general  linear  programming  problem  of  the  form:  Find 

x,  ,x„,x.,. . . ,x  ,  such  that 
12  3  n 


n 


j-1 


(1) 


is  a  minimum  over  all  sets  x, , x_ , x„ , . . . , x  that  satisfy 

12  3  n 

the  independent  constraints: 

n 

I  aijXj  =  bi  1=1 , 2 , 3 , . . . ,rn  ,  (2) 

j-1 


and  x.  *  0  j=l,2,...,n,  where  C.,a..,b.  are  fixed  numbers. 

J  J  ij  i 


USAGE 


a 

C 


ij 


The  a.  .  is  stored  in 
ij 

in  cell  A(i, j ) ;  is 
with  Cj  in  cell  C(j); 


a  two-dimensional  array,  A,  with 
stored  in  a  one- d imens ional  array, 
and  b^  is  stored  in  a  one-dimension 
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The  calling  sequence  is 

i  LL  SIMPLE  (II,M,N,A,B,C,KO,X,P,JH,XX,Y  PE,E) 

where 

II  =  0; 

M  =  Number  of  rows,  m; 

N  =  Number  of  variables,  n; 

A,  B,  C  are  as  above; 

KO  =  A  subscripted  variable 
of  dimension  6; 

X  =  A  subscripted  variable 
of  dimension  n  or  more; 

P,  JH,  XX,  Y,  and  PE  =  Subscripted  variables  of 

dimension  m  or  move,  and 

E  =  A  subscripted  variable 
of  dimension  m^  or  more. 

The  dimension  of  A  (line  0008  of  FORTRAN  listing)  muse 
agree  (at  least  in  the  first  subscript)  with  the  dimension 
of  A  in  the  calling  program.  The  other  dimensions  need  not 
agree  with  those  of  the  calling  program.  The  subroutine 
does  not  change  II,  M,  N,  A,  B,  or  C. 

RESULTS 

Upon  exiting  from  the  subroutine, 

X(l) ,X(2) , . . . ,X(n)  contains  x,,x0,...,x  (the  solution); 

i  z  n 

P(l) , P(2) , . , , , P(m)  cont-ins  the  shadow  prices,  i.e., 

the  negative  of  th _•  dual  solution; 
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K0(1)  contains  an  0  if  the  problem  was 
feasible,  1  if  the  problem  was 
infeasible,  2  if  the  problem  had 
an  infinite  solution,  and  4  or  5 
if  the  algorithm  did  not  terminate; 

K0(2)  is  the  number  of  iterations  taken; 

K0(3)  is  the  number  of  pivots  performed 
since  the  last  inversion; 

K0(4)  is  the  number  of  inversions  per¬ 
formed  ; 

K0(5)  is  the  number  of  pivot  steps  per¬ 
formed  ; 

and 

K0(6)  contains,  if  the  problem  has  an 
infinite  solution,  the  number  of 
the  variable  that  was  infinite. 


If  an  initial  basis  is  available,  this  basis  may  be 
communicated  to  the  subroutine  by  letting 

II  =  1  , 

I  0.0  if  variable  i  is  not  in  basis, 

X(i)  «< 

I  (non-zero)  if  variable  i  is  in  basis, 

and  the  other  quantities  remain  as  above. 

If  the  constraints  (2)  are  linearly  dependent,  the 
problem  is  considered  infeasible.  II  the  problem  is  in¬ 
feasible,  (K0(1)  =  1),  the  P(i)  contains  the  coefficients 


-15- 


that  generate  the  infeasible  row;  i.e.,  if  the  problem 
is  infeasible,  the  infeasibility  is  represented  by 


n  m 

£  D( j )X( 1  P(i)B(l)  =  E 

j-1 


where 


m 

D(j)  =  £  P(t)A(t,j)  j-l,2,...,n  . 

i  =  l 


Then,  when  infeasible,  D( 1) , D(2) , . . . ,D(n) ,  and  -Ewill 
simultaneously  all  be  either  non-negative  or  non-positive; 
and  if  E  is  zero,  then  D(l) ,D(2) , . . . ,D(n)  are  all  zero. 

If  the  problem  has  an  infinite  solution,  (K0(1)  =  2), 
define  a  vector  XI  as  follows: 

D0  1  J=1,N 

1  X1(J)  =  0. 

D0  2  1=1, M 

J=  JH( I) 

2  X1(J)  =  -Y (J ) 

J  =  KO  (6) 

X1(J)  =  1. 


Then  the  infinite  solution  is  represented  by 


lim  (X(J)  +  t-Xl(J)) 

t-® 
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If  the  problem  did  not  terminate,  then  4m+10  itera¬ 
tions  were  insufficient  to  complete  the  solution.  If 
K0(1)  =  5,  the  problem  is  not  yet  feasible;  and  if  K0(1)  =  4, 
the  problem  is  feasible  but  not  yet  optimal.  The  subroutine 
may  be  called  a  second  time,  with  11=1,  or  the  equation  on 
card  0021  may  be  changed  to  allow  for  more  iterations. 

EXAMPLE 

Let 


3xl 

+ 

14x.  +  lxc 
4  5 

+ 

II 

X 

7  , 

lx2  + 

16x^  +  ^x5 

- 

2x6  ■ 

5  , 

lx^  + 

lx4 

= 

0  , 

and 

minimize : 

-  28x4  - 

x5 

-  2x, . 
0 

Using  the  subroutine  with  DIMENSION  A(50,136),  we  zero 
out  A  and  set 


A(l.l) 

= 

3. 

C(l) 

=  0. 

A(2,2) 

= 

1. 

C(2) 

=  0. 

A(3 , 3) 

= 

1. 

C(3) 

=  0. 

A(l,4) 

« 

14. 

C(4) 

=  -28  . 

A(2,4) 

= 

16. 

A(3,4) 

= 

1. 

A(l,5) 

= 

1. 

C(5) 

=  -1  . 

A(2,5) 

= 

0.5 

A(3.5) 

= 

1. 

A(1 , 6) 

ss 

1. 

C(6) 

=  -2  . 

A(2,6) 

ss 

-2. 

B  ( 1) 

= 

7. 

B(2) 

= 

5. 

B(3) 

= 

0. 

M 

= 

3 

N 

= 

6 

II 

- 

0 
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Then  we 


CALL  SIMPLE  (II ,M , N , A , B ,C , KO ,X , P , JH ,XX , Y , PE, E) 


where  B,P,JH,XX,Y,  and  PE  have  DIMENSION  of  at  least  3 
(i.e.,  M) ,  KO  has  DIMENSION  6,  C  and  X  have  DIMENSION  of 
at  least  6  (i.e.,  N) ,  and  E  has  DIMENSION  of  at  least  9 
(i.e. ,  M2) . 

Upon  returning  from  the  subroutine,  we  find  the  fol¬ 
lowing  values  in  storage: 


K0(1)  =  0 
KO  (2)  =  2 
KO(3)  =  2 
KO  (A)  =  1 
KO(5)  =  5 
K0(6)  =  0 


Optimum  solution, 

Number  of  iterations  taken, 

Number  of  pivots  since  last  inversion, 
Number  of  inversions  performed, 

Total  number  of  pivots, 

Not  applicable. 


Also 

P(l)  -  2. 

P(2)  =  0.  , 

P(3)  ■=  2.38  x  10  0. 


and 

X( 1)  =  o. 

X(2)  -  18.9999998  19. 

X(3)  -  0. 

X(4)  -  -0 . ■  0. 

X(5)  -  0. 

X(6)  -  6.99999994  s  7. 
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Now  we  may  calculate  the  reduced  costs  from  this  informa¬ 
tion  using,  for  example,  the  following  instructions  to  put 
t  h 

the  j —  reduced  cost  into  CBAR(j): 

DO  12  J  =  1,N 
CBAR(J)  =  C(J) 

DO  11  I  =  1 ,M 

11  CBAR(J)  =  CBAR(J)  -I-  P(I)*A(I,J) 

12  CONTINUE 

This  results  in 

CBAR(l)  =  6. 

CBAR(2)  =0.  7 

CBAR(3)  =  2.384  x  10" /  -  0. 

CBAR(4)  =  2.384  x  10~7  a  0. 

CBAR(5)  =  1. 

CBAR (6)  =  0. 

Following  is  a  listing  of  the  FORTRAN  program  written 
in  FORTRAN  IV  [4 ',  then  an  ALGOL  procedure  using  McCracken 
[5].  The  ALGOL  procedure  is  a  translation  of  the  FORTRAN 
program,  so  any  discrepancy  should  be  resolved  in  favor  of 
the  FORTRAN  program.  Table  1  gives  the  names  and  meanings 
of  the  symbols  used  in  the  calling  sequence  to  the  FORTRAN 
pregram,  ir  the  FORTRAN  program  itself,  and  in  the  ALGOL 
procedure . 


SYMBOL  NAMES  AND  MEANINGS 


* 

* 
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"O 

C  • 
<0  C 
C 
N 

w  3 


o  o 

Xi  cn 
G 

>.  G 
U)  si 
u 


X  -I 
O  w 
CJ  o 
X  &c 

<C  t— i 

oo 


QJ 

SI  o 


u  x 


o 

iJ  U-4 

o 

oo 

C  (A 

•r4  01 

TO  60 
C  (0 
O  4-1 
Q.  W 

(0 

04  4J 

Vj  C 
V-i  04 
O  J-i 
CJ  01 
u-i 
J)  U-i 
04  -r4 
—I  XJ 


•r-4  4-1 
4J  CO 
C 

C0  u 
3  3 
cr  xi 

04  - 

X  c 

4-1  9 

•r^ 

'  u 
04  3 
C  CJ 
v-4  9 

U  r— 4 

3 

9  04 

V4  G 
X  co 
3  W 
CA 

04 

^0  X 

H  C 

QS 

O 

U*  X3 
04 

0)  V-i 

X  9 

4J  U 
Cfl 
C 

M  04 
CO  u 
CO 


cq 
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FORTRAN  PROGRAM 


•ibftc  simple  rtf 

C  AUTOMATIC  SIMPLEX  REDUNDANT  EQUATIONS  CAUSE  INFEASIBIL l T V 

SUBROUTINE  SIMPLE  I  I NFl AG .N* .NN .A  ,  B . C ,RO .KB ,P . JM , X . V .PE . E I 

real  bu  i  .c.  1 1  ,pi  1 1  .x<  1 1  ,rm  .pei  1 1  .e  <  1 1 

INTEGER  I  NFL  AG. NI.NN.RO I  6 ) • KB ( 1 1 » JH  I 1) 

EQUIVALENCE  IXX.LLI 

C  THE  FOLLOWING  DIMENSION  SHOULD  BE  THE  SAME  HERE  AS  IT  IS  IN  CALLER. 
REAL  A ( SO . 1  SB  I 

real  AA.AI JT.BB.COST.DT.RCOST.TExP.TPIV.TY.XOLD.XX.XV.rl ,YMAX 
INTEGER  I . I  A . INVC • IR. 1  TER. J. JT .K.XB J.L .LL •M.M2.MM.N 
INTEGER  NCUT . np I  V .NUMVR .nve r 
LOGIC a;  FEAS.VER.NEG.TRIG.RQ. absc 

c 

C  SET  INITIAL  values.  SET  CONSTANT  VALUES 

ITER  ■  0 
NUMVR  ■  0 
NUMPV  •  0 
M  ■  MX 
N  ■  NN 

TEXP  •  . S#*  IB 

NCUT  •  4»M  ♦  10 

NVER  •  M/2  ♦  S 

M2  ■  M«*2 
FEAS  •  .FALSE. 

IF  ( INFLAG. NE. 01  GO  TO  1400 

C*  'NEW  START  phase  one  with  singleton  basis 
DO  1402  J  ■  l.N 
RBIJI  ■  0 
RQ  •  .FALSE. 

DO  1401  I  •  l.M 

IF  I  At  I . J) .EG. 0.01  GO  TO  1401 
IF  IKQ.OR.AI  I.Jl.tT.O.Ol  GO  TO  1402 
RQ  •  .TRUE. 

1401  CONTINUE 
RBIJI  ■  1 

1402  CONTINUE 

1400  DO  1401  I  •  1  «M 

JH  (I)  •  -1 

.4Q1  CONTINUE 

C*  •  VER '  CREATE  INVERSE  FROM  ' KB •  AND  *JH*  (STEP  71 

1120  VER  •  .TRUE. 

INVC  ■  0 

NUMVR  ■  NUMVR  *1 
TRIG  ■  .FALSE. 

00  1101  |  «  1.M2 

Fill  ■  0.0 
1101  CONTINUE 
MM.  1 

DC  1 1 11  I  •  1  .M 
El  MM)  ■  1.0 
PE  I  I  I  ■  0.0 

XI I)  •  Bill 

IE  I  JH  I  I  I  .NE.OI  JHII)  ■  -1 

MM  •  MM  ♦  M  ♦  1 

11 1 1  CONTINUE 


L  SUB0001 
LSUB0002 
L  SUB0001 
LSUB0004 
LSUB0005 
LSUB0006 

lsubooot 

LSUB0008 
LSUB0009 
LSUBOOl 0 
LSUB001 1 
L  SUBOO 1 2 
LSUBOOl 3 
L  SUBOO 14 
L  SUBOO 1 S 
L  SUBOO 16 
L  SUBOO 1 7 
LSUBOOl 8 
L  SUBOO 1 9 
LSUB0020 
L  SUBOO? 1 
LSUB0022 
LSUB0023 
LSUB0024 
LSUB002S 
L SUBOO 26 
L  SUB002  7 
L  SUB0028 
LSUB0029 
LSUB0030 
L  SUBOO  3 1 
L  SUBO0 12 
LSUB0033 
L  SUBOO 14 
LSUBE/IS 
L  SUBOO  3fc 
L  SUBOO  3  7 
L  SUBOO  38 
L  SUBOO  39 
LSUB0040 
L  SUB0041 
LSUB0042 
LSUB0043 
L  SUBOO 4* 
L  SUB004S 
LSUB0046 
LSUB0047 
LSUB0048 
L  SUB0049 
LSUB0010 
L  SUB0011 
LSUBOOl? 
LSUB0013 
L  SUB0014 
L  SUBOOSl 
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c 

FORM  INVERSE 

L  SUB00S6 

00  110?  JT  ■  l.N 

L  SUB005  7 

IF  (K6I JTI.EO.OI  GO  TO  U5? 

lsuboosb 

GO  TO  600 

L  SU600S* 

C  600 

CALL  JMY 

L  SUB006C 

c 

CHOOSE  PIVOT 

L  SUB0061 

111* 

TY  ■  0.0 

LSUBOOo? 

KO  ■  .FALSE. 

LSUB006* 

oo  no*  i  ■  i  (M 

LSUB006* 

IF  I  JHI  1  1.NE.-1.0R.ABSI  Yl  IM.LE.TPIVI 

GO  TO 

no* 

LSUB006* 

IF  (KOI  GO  TO  1116 

LSUB0066 

IF  (XI I I.EO.O.I  GO  TO  111* 

L  SUB006T 

IF  IABSI Yt  1  l/Xt 1 1 I.LE.TYl  GO  TO  110* 

L  SUB0066 

TY  •  ABS  (  Y  (  Il/Kll  1 

LSUB006* 

GO  TO  11  IB 

LSUB0070 

1m 

KO  >  .TRUE. 

lsubooti 

GO  TO  1117 

LSUB007? 

1116 

IF  I  X  I 1 1 .NE.0..OR.ABSIY 1 1 1 }  .LC • T Y 1  GO 

to  no* 

lsubooti 

1117 

TY  •  ABS  <  Y ( 1  I  1 

lsuboot* 

111* 

IR  •  1 

LSUB007* 

11C* 

CONTINUE 

L  SUB0076 

KB  I JT I  •  0 

LSUB0077 

C 

TEST  PIVOT 

LSUB007B 

IF  1  TY.LE.O. 1  GO  TO  110? 

lsuboot* 

c 

PIVOT 

lsuboobo 

GO  TO  *00 

lsuboobi 

C  *00 

call  piv 

LSUBOOB? 

110? 

CONTINUE 

LSUBOOBI 

c 

RESET  ARTIFICIALS 

LSUBOOB* 

oo  no*  i  •  i.n 

LSUBOOB* 

IF  I JH( 1 i.EO.-l  1  JHMl  ■  0 

L  SUB00B6 

IF  «JH»n.E0.01  FEAS  ■  .false. 

LSUBOOB 7 

1  10* 

CONTINUE 

LSUBOOBB 

1?00 

VER  •  .FALSE. 

LSUBOOB* 

e 

•••  PERFORM  one  ITERATION 

••• 

LSUB00*0 

C*  'XCK*  DE  T  E RN I NE  FEASIBIL  1  TY 

(STEP 

1 1 

L  SUB00*1 

NEG  •  .FALSE. 

LSUBOO*? 

IF  IFEASI  &0  TO  *00 

LSUBOO** 

FEAS-  .TRUE. 

LSUBOO** 

00  1201  |  •  l.N 

LSUBOO** 

IF  IXIII.LT. 0.0)  00  TO  1?*0 

LSUB0096 

IF  (JHlII.EO.Ot  TEAS  •  .FALSE. 

L  SUB009T 

1201 

CONTINUE 

lsuboo*b 

C*  »GET  •  GET  APPLICABLE  PRICES 

1  STEP 

?i 

LSU6009* 

IF  I.NOT.FEAS)  GO  TO  *01 

LSUB0100 

*00 

00  *01  1  •  l.N 

LSUB0101 

PI  II  •  PEI  1  1 

LSUB010? 

IF  IXIII.LT. 0.1  II  II  •  0. 

LSUB010) 

*01 

CONTINUE 

LSUB010* 

ABSC  ■  .FALSE. 

LSU6010* 

GO  TO  *** 

L  SU0O1 06 

1?*0 

FEAS  •  .FALSE. 

LSU60107 

NEG  •  .TRUE. 

LSUB010S 

*01 

00  *0*  J  ■  1  •  N 

LSUB010* 

PIJI  ■  0. 

LSUB0110 
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•>0* 

CON  T  I NUE 

ISUB0111 

ABSC  •  .TRUE. 

LSUB0112 

DO  505  I  •  1  .M 

LSUB0113 

MM  •  1 

LSUB0U4 

IE  IX! 1 1 .GE.0.01  GO  TO  507 

LSUB0113 

ABSC  •  .FALSE. 

LSUB01 16 

DO  500  J  •  l.M 

L  SUB01 1 7 

PIJI  -  PIJI  ♦  E 1  MM ) 

L  SUB01 10 

MM  •  MM  ♦  M 

LSUB0119 

500 

CONTINUE 

L  SUB01 20 

GO  TO  505 

LSUB0121 

507 

IF  1  JHI I I.NE.Ol  GO  TO  505 

LSUB0122 

IF  (Kill. NE. 0.1  ABSC  •  .FALSE. 

LSUB0123 

DO  510  J  •  1  .M 

LSUB0124 

P  1  Jl  •  PIJI  -  E 1  MM ) 

L  SUB0123 

MM  ■  MM  ♦  M 

LSUBO120 

510 

CONTINUE 

L  SUB01 27 

505 

CONTINUE 

LSU6O120 

C •  ‘MIN'  FINO  MINIMUM  REOUCED  COST 

(STEP 

3) 

LSUB0129 

509 

JT  ■  0 

LSUB0130 

BB  •  0.0 

LSUB0131 

DO  701  J  ■ 1 » N 

LSUB0132 

IE  1  KB  I J 1 . NE.O 1  GO  TO  701 

LSUB0133 

DT  •  0.0 

LSUB0134 

DO  303  I  •  l.M 

L  SUB01 35 

DT  ■  DT  ♦  Pill  •  A  1 1 • J 1 

LSUB0136 

30  3 

CONTINUE 

L  SUB01 37 

IF  (FEASI  DT  •  DT  ♦  CIJI 

LSUBO130 

IF  1 ABSC  1  DT  •  -  ABSIOTI 

LSUB0139 

IF  (DT.GE.BBI  GO  TO  701 

LSUB0140 

BB  •  DT 

LSUB0141 

JT  •  J 

LSUB0142 

TO  1 

CONTINUE 

LSUB0143 

C  TEST  FOR  MO  PIVOT  COLUMN 

LSUB0144 

IF  (JT.LE.0I  GO  TO  203 

LSUB0145 

C  TEST  FOP  ITERATION  UNIT  EXCEEDED 

LSUB0146 

IF  I  1 TER.GE.NCUT |  GO  TO  160 

LSUB0147 

ITER  ■  ITER  ♦  1 

LSUBO140 

C  •  •  JMY •  MULTIPLY  INVERSE  TIMES  AI..JT) 

(STFP 

*  1 

L  SUB01 49 

600 

DO  610  1 •  l.M 

LSUB0150 

Y( 1 1  •  0.0 

L  SUB0131 

610 

CONTINUE 

LSUB015? 

LL  •  0 

LSUB0153 

COST  -  CIJTI 

LSUB0154 

DO  605  1 «  l.M 

L  SUB01 55 

Al JT  ■  All*  JT  > 

LSUBO150 

IF  ( Al JT.EQ.O.  1  GO  TO  60? 

LSUB0157 

COST  •  COST  ♦  A 1 JT  •  PE  (  1  I 

L  SUB01 50 

DO  606  J  •  l.M 

LSUB0159 

t  L  ■  L  L  ♦  1 

L  SUB0160 

YIJ)  ■  Y 1 J 1  ♦  Al JT  •  EILLI 

LSUB0161 

606 

CONTINUE 

LSUBO10? 

GO  TO  605 

LSUB0163 

60? 

LL  •  LI  ♦  M 

L  SUBO104 

603 

CONTINUf 

L  SUB0165 
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C  COMPUTE  PIVOT  TOLERANCE 

YMAX  »  0.0 
00  6? 0  I  *  1  .M 

YMAX  ■  AMAXH  ABS  (  Y  (  I  *1  .  YMA  X  ) 

620  CONTINUE 

TP  t  V  ■  YMAX  •  TExP 

C  RETURN  TO  INVERSION  ROUTINE,  IE  INVERTING 

IE  (VERl  GO  TO  111* 

C  COST  TOLERANCE  CONTROL 

RCOST  ■  YMAX/BB 

IE  (TRIG. AND. BB.Gf.-TPIVI  GO  TO  20* 

TRIG  •  .EALSE. 

IE  IBB.GE.-TPIv!  TRIG  •  .TRUE. 

C*  'ROW •  SELECT  PIVOT  ROM 
C  AMONG  EQS.  WITH  X  *0 •  F  I  NO  MAXIMUM  Y 

C 


(STEP  S) 

AMONG  ARTIFICIALS.  OR. 


IE  NONE. 


LSUB0166 
LSUB0167 
LSUB016B 
L SUBO 169 
LSUB01 7 0 
LSUB01 71 
L  SUB01 7? 
L  SUB01  71 
LSUB01 74 
LSUB01 7S 
L  SUB01 76 
LSUB01 77 
LSUB01 78 
L  SUB01 79 
LSUB01 80 


GET 

MAX 

POSITIVE  Y ( I »  AMONG  REALS. 

LSUB0181 

IR  ■ 

0 

LSUB0182 

AA  ■ 

0.0 

L  SUBO 18  3 

KO  - 

. F  AL  SE • 

L  SUB01 84 

no 

1030  I  ■  1  ,  M 

L  SUBO 183 

IF 

I  X(  I  1  .NE.C .O.OR. Y « 1 l .LE.TPI V 1  GO  TO 

1030 

LSUB01 86 

IF 

( JH ( I l.EO.Ol  GO  TO  1044 

L  SUBO 187 

IF 

IKQI  GO  TO  1030 

L  SUBO 1 88 

1043 

IF 

(Y(  1  t.LE.AAl  GO  TO  1030 

L SUBO 189 

GO 

TO  1047 

L  SUBO 190 

1044 

IF 

IKQI  GO  TO  1043 

L  SUBO 1 9 1 

KO 

•  .TRUE. 

L SUBO 192 

1047 

AA 

■  V  f  1  1 

L  SUBO 193 

IR 

•  1 

L SUBO 194 

1030 

CONTINUE 

LSUB0193 

IF  ( IR.NE.OI  GO  TO  1099 

L SUBO 196 

AA  • 

1 .0E»20 

LSUB0197 

FIND  MIN.  PIVOT  AMONG  POSITIVE 

equations 

L  SUB01 98 

00  1010  1  •  l.M 

L SUBO 199 

IF 

(Y(I1.LE.TPIV.OR.K(||.LE.O.O.OR.Y(I)*AA.LE.K(I> 

»  GO  TO  1010  LSUB0200 

AA 

•  iiii/yiii 

LSU80201 

IR 

■  i 

LSU60202 

1010  CONTINUE 

IF  ( . NOT . NEG I  GO  TO  1099 

C  F I  NO  PIVOT  AMONG  NEGATIVE  EQUATIONS.  IN  WHICH  X/Y  IS  LESS  THAN  THE 
C  MINIMUM  X/Y  IN  THE  POSITIVE  EQUATIONS.  THAT  HAS  THE  LARGEST  A6SF ( Y I 
BB  •  -  TP  I V 
DO  1030  I  ■  l  , M 

IF  ( X(  1  I  .GE .0..OR.Y l  I  I .GE .8B.0R.YI I « 4AA.GT .X  I  I  I  I  GO  TO  1030 
BB  ■  Y( I l 
I R  •  I 
1030  CONTINUE 
C  TEST  FOR  NO  PIVOT  ROW 
1099  IF  ( IR.LE.Ol  GO  TO  207 
C •  • P 1 V  •  PIVOT  ON  (  |R, JT  » 

I  A  •  JH(  | R I 

IF  < I A.GT.OI  XB ( IA I  -  0 
900  NUMPV  •  NUMPv  ♦  1 

JH ( I R I  •  JT 
KB ( JT  )  •  | R 


(STEP  6) 


L  SUB0203 
LSUB0204 
LSUB0203 
L  SUB0206 
LSUB0207 
L  SUB0208 
LSUB0209 
LSUB0210 
LSUB021 1 
L  SUBO? 1 2 
LSUB0213 
LSUB0214 
LSUB0213 
LSUB0216 
L  SUB0  2  1  7 
LSUB0218 
LSUB02I9 
L  SUB0270 
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VI  •  -Y( IR  t 

L  SUB022  1 

VI IR»  •  -1.0 

L  SUB0222 

LL  ■  0 

L  SU0O22  3 

transform 

INVERSE 

LSUB0224 

DO  904  J  •  l.M 

LSU80225 

L  •  LL  ♦  IR 

L  SU60226 

IF  1 E 1 L 1 .RE .0.0 1  GO  TO  90S 

LSUB0227 

LL  •  LL  ♦  H 

LSUBO220 

GO  TO  904 

LSUB0229 

905 

XY  ■  EIL  I  /  VI 

LSUB023C 

PEI J1  •  PE  1 J»  ♦  COST  •  XY 

LSUB0231 

EIL  1  ■  0.0 

L  SUB02  3  2 

DO  906  l  ■  1  .M 

LSUB0233 

LL  •  LL  ♦  1 

LSU60234 

EILLI  ■  EILLI  ♦  XV  •  V< || 

LSUB0215 

906 

CONTINUE 

LSUBO230 

904 

CONTINUE 

LSUB0217 

TRANSFORM 

X 

LSUBO230 

XV  ■  X 1 1 R 1  /  VI 

LSUB0239 

DC  900  I  -l.M 

LSU60J40 

XOLD  •  X  1  1  1 

LSUB0241 

XIII  •  XOL0  ♦  XV  •  V  (  1  1 

LSUB0242 

IF  I.NOT.VER.AND.XI I t.LT.O.*  AND . 

XOLD.GE.O. 1  XI  1  I  ■  0. 

LSUB0243 

900 

CONTINUE 

L  SUB0244 

VI  IRl  •  -VI 

LSUB0245 

X 1 I R |  «  -XV 

LSUBO240 

IF  IVERl  GO  TO  1102 

LSUB0247 

IF  INUMPv.LE.Mi  GO  TO  1200 

LSUBO240 

TEST  F OR  INVERSION  OR  THIS  ITERATION 

L  SUB0249 

I N VC  •  I NVC  41 

LSUB0250 

IF  I INVC.EO.NVERI  GO  TO  1320 

LSUB0251 

GO  TO  1200 

LSU60252 

END  OF  ALGORITHM.  SET  EXIT  VALUES 

•  •• 

LSUB0253 

207 

IF  ( .NOT.Ff AS.OR.RCOST.LE.-IOOO. ) 

GO  TO  2u  3 

LSUB0254 

INFINITE  SOLUTION 

LSUB0255 

X  -  2 

L  SUBO250 

GO  TO  250 

LSUB0257 

PROBLEM  IS  CYCLING 

LSUBO250 

160 

X  ■  4 

L  SUB0259 

GO  TO  250 

LSUB0260 

FEASIBLE  OR  infeasible  SOLUTION 

L  SUBO20 1 

203 

X  •  0 

LSU8O202 

250 

|F  I.NOT.FEASI  X  ■  X  ♦  1 

LSUB0263 

DO  1399  J  •  l.N 

LSUBO204 

XX  •  0.0 

LSUB0265 

XBJ  .  XBIJI 

LSUBO200 

IF  IXBJ.NE.O)  XX  ■  XUBJI 

LSUBO207 

XBIJI  •  LL 

LSUBO260 

399 

CONTINUE 

L  SUBO209 

XOI 1 1  ■  X 

LSUB02T0 

XOI2I  •  ITER 

L  SUB02  7 1 

XOI3I  ■  INVC 

LSU602T? 

XO 1 4 )  •  NUMVR 

LSUB02T3 

XOI 51  ■  NUMPV 

LSUB0274 

XOI6I  ■  JT 

LSUB02T5 

RETURN 

LSUB0276 

END 

LSUB0277 

277 

277 
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ALGOL  PROCEDURE 
Ok  used  in  place  of  X ) 


procedure  SIMPLF  ( in  f  lag,  ir ,  n,a  ,  b ,  c ,  ko ,  kb ,  p ,  jh ,  x ,  y ,  f  e ,  e ,  z)  ; 
value  m,  n;  rea  1  i  rra  y  hf  p,  x ,  y  #pe  ( 1 :  nr)  ,c ,  z  [  1  :  n)  ,  e  ( 1 :  mt  2)  , 
a  ( 1 :  m,  1 :  n)  ;  i  nteger  ar  ray  ko  f  1 :  b)  ,  jh  [  1 :  ir |  ,  k  h  ( 1 :  n)  ; 
integer  m,n, inflaq; 

coirment  SIMPLE  solves  the  linear  programming  protlem: 

find  z(1)  #  z  (2)  ,  z  (3)  , . . .  ,  z  (n)  that 
n 

minimizes  l  c(j)  •  z  (j)  subject  to 

j=1 

n 

£  a(i,j)  *  z  ( j)  =  b  (i)  for  i=  1,  2,  3, . . .  f  ro  and 
1=1 

x  (1)  >0  #  x  (2)  >0,  x(3)>0,...,  x  (n)  ^0. 

SIMPLF  simultaneously  solves  the  dual  problem: 
find  w  (1)  ,  w  (2)  ,  w  (  3)  ,  . . .  ,  w  (m)  that 
m 

maximizes  l  b  (i)  •  w  (i)  subject  to 

i=  1 

m 

£  a  (i ,  j)  ♦  w  (i)  <  c(j)  for  j=  1  ,  2,  3, . . .  ,n 

i  =  1 

wh«  re  no  sign  restriction  is  put  on  w.  The  NEGATIVE 
of  w  appears  in  p. 

inflag,  m,  n,  a,  b,  and  c  are  input  quantities:  they  are 
not  chanqed  by  the  procedure.  The  meaninqs  of  it,  n,  a,  fc 
and  c  are  as  given  in  the  above  equations.  inflaq  is  nor 
mally  input  as  a  zero:  this  signals  the  procedure  that  no 
initial  l*asis  is  provided.  If  inflaq  is  input  as  a  non¬ 
zero  number,  the  procedure  will  expect  a  basis  to  be  pro¬ 
vided  by  setting: 

kb(j)  =0  if  column  j  is  not  in  basis, 
and  kb(j)  #0  if  column  j  is  in  the  basis. 
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Atter  the  ; rocedure  is  finished,  ko  (1)  contains  a  number 
that  denotes  the  condition  of  the  problem,  and  z  and  p 
contain  the  numeric  answers. 

If  the  problem  is  feasible  and  optimal,  then  ko  (1)  =  0. 

If  the  problem  admits  no  feasible  solution,  then  ko(1)  =1, 
and  if  the  problem  has  an  infinite  solution  (i.e.  dual 
infeasible),  then  ko  (1)  =  2.  If  the  alqorithm  does  not 

terminate  after  4m*10  iterations,  ko  ( 1 )  is  set  to  4  if  the 
problem  is  feasible,  and  to  5  if  yet  no  feasible  solution 
has  been  found; 
begin 


real  ^  a,bb,  cost ,dt, rcost, texp, tpi v, ty,ymax; 
integer  i,invc,ir,iter,  j,jt,k,m2,  rrm,  ncut ,  npiv,numvr , 
nver;  boolean  absc, f eas , kq, neq, tr ig , ver , 
iter:=  6; 
numvr:=  0; 
numpv:=  0; 
texp: =  2. t  (- 16) ; 
ncut:3  4*m  ♦  10* 
nver  :=  m/2  ♦  S; 

m2:=  mt2; 
feas:=  false; 
if  (inflaq  *  0)  then 


comment  start  phase 

for  j : =  1  step  1  until 
n04 :  begin 

STFTnr :  =  0; 

kq: -  false; 

for  i:=  1  step  1  until  m  do 
n05:  begin 

if  (a  (i  ,  j)  =  0.)  then 
i_f  (kq  V  a  (i ,  j)  <  0.) 


go  to  nOO; 
one  wFth  singleton  basis; 
n  do 


go  to  n 0  3 ; 
then  go  to 


n02 


n03 : 


kq : =  true ; 
end  nOS; 
kb  ( j)  :=  1; 


n02:  end  n04  ; 

nOO:  for  1:=  1  step  1  until  m  do  jh  (i)  : =  - 1 ; 
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cororent 
m20 : 


create  inverse  from  "kb"  and  *  jh" 


(step  7)  ; 


k20 : 


ver:*  true; 
invc:=  0; 

mimvr  :=  numvr  ♦  1 ; 
trig;=  false; 

for  i ;  =  1  step  1  until  m2  do  e  ( i )  :  =  0.0; 
mm:  =  1 ; 

for  i : =  1  step  1  until  m  do 
begin 

e  (mm)  :  =  1.0;  pe(i):=  0;  x(i):=  b  (i)  ; 
if  (jh(i)  *  0)  then  jh(i):=  -1; 
mm:3  mm  ♦  m  ♦  1; 


k  1  3 :  end  k  2  0 ; 

for  jt : =  1  step  1  until  n  do 
k  2 1 :  begin 

if  (kb(jt)  =  0)  then  go  to  k02; 

Get  Column  ( jt, a, c,  e, pe, y ,m, cost , texp, tpiv ,  ymax)  ; 
ty:  =  0.0;  kq:=  false; 
for  i : =  1  step  1  until  n  do 
k22:  begin 

if  ( jh  (i)  *  -1  V  abs  (y  (i)  )  <  tpiv)  then  go  to  kOk 

if  (kg)  then  go  to  k16; 

if  (x  (i)  =  0)  tTTen  go  to  k15; 

O  (ahs  (y  (i)  /x  (i)  )  <  ty)  then  go  to  k0*; 

ty:  =  abs  (y  (i)  /x  (i)  )  ; 
go  to  k  1  8  ; 
k15:  kq:-  true; 

<12  Ifl  k  1  7  ; 

k  1 6 :  JLf  (x  (i)  i  0.0  V  abs  (y  (i)  )  >  ty)  then  go  to  kOk; 

k  1  7:  ty :  =  abs  (y  (i)  )  ; 

k18:  ir:  -  i; 

k04:  end  k22; 

kb(jt)  :=  0; 

if  (ty  <  0.0)  then  go  to  k02; 

Pivot  (ir , jt ,e, jh, kb, pe,x,y#m,cost, numpv , ver)  ; 
k02:  end  k21; 

for  i : =  1  step  1  until  m  do 
begin 

if  (jh(i)  *  -1)  then  jh  (i)  :*  0; 

if  ( jh  (i)  =0)  then  feas:*  false 

end; 

100:  ver:*  false; 

comment  determine  feasibility  (step  1) ; 

neg : *  false ; 

if  (feas)  then  go  to  eOO; 

Tea s :  true ; 

for  i:=T~step  1  unti  1  rc  do  _if  (x  (i)  <0.)  then  go  to  150 
else  if  (jh  (1)  3  then  feast*  false  ; 
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comroent  get  applicable  prices 

if  (  ifeas)  then  30  to  eO  1 ; 
eUO:  for  i:=  1  step  1  until  ir  do 

~~begi  n 

p  (i)  :=  [H*  (i)  ; 

if  (x  ( i)  <  0.0)  then  x  (i)  :  3  0.0 

end ; 

a  bs  c :  =  false ; 
go  to  eT9 ; 

150:  feas:=  fa lse; 

nog  :=  true; 

e0 1 :  for  j:  3  1  step  1  unt i  1  rr.  _do  p  (j)  :  =  0.0 
absc:3  true; 

for  i :  3  T  step  1  unt i  1  m  dio ; 
ell:  beg i n 

mm: =  i; 

if  ( x  ( i )  >  0.0)  then  go  to  e07; 

absc: 3  false; 

for  j  :  =  1  step  1  unt i 1  m  do 
becjin 

P  (1)  :=  P  ( j)  ♦  e  (irm)  ; 
mm :  =  mm  ♦  m 
end ; 

go  to  e05 ; 

e07 :  if  TTh  (i)  *  0)  then  <jo  to  e05; 

IT  (x (i)  #0.)  then  absc:3  false; 

for  j:=  1  step  1  until  n  do 
bec[  i  n 

P  (1)  •*=  F  ( j)  -  e  (mm)  ; 
mm :  3  mm  ♦  m 
end ; 

e05 :  end  e  11; 

comment  find  minimum  reduced  co3t 

e9  V:  jt :  =  0; 

bb  :=  0.0; 

g0  2:  for  j:=  1  step  1  until  n  do 
begin 

if  (kb(j)  #  0)  then  go  to  g01; 

It :  -  0.0; 

for  i :  =  1  step  1  until  m  do  dt :  =  dt 
if  (feas)  then  dt:3  dt  ♦  c  ( j)  ; 
if  (absc)  then  dt:=  -abs  (dt)  ; 
if  (dt  ?.  bb)  then  go  to  g()1  el Sf  bb: 
jt:*  l; 

g0  1  :  end  g02  ; 

if  ("jt  s  0)  then  qo  to  bOI; 

Tf  (iter  >  ncut)  then  30  to  a  60; 
iter:=  iter  ♦!; 


(step  2) 


(step  3) 

P  (i)  (i »  j)  ; 
d  t ; 
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cotr merit  multiply  mvorsn  times  d  (.  ,  )t)  (stef  4)  ; 


fOO: 

Get  Co  lumn  ( jt  ,d  ,  c,  e  ,  pe#  /#mrcost,texp,tpiv, 
roost :  =  ynid x/H  ; 

if  (trig  A  hi  -*  -tr  iv)  then  go  to  HH; 
trig:-  f  d  lse; 

if  (hb  >  -tpiv)  then  trig:=  true; 

yndx)  ; 

comment  select  fivor  row 

i  r  :  =  0  ; 

d  u  :  —  0.0; 

kg: =  f d lse; 

for  i:=  1  stef  1  until  v  io 

(Stef 

102  : 

be^g  i  n 

if  (x  (l)  #  0.0  V  y  (i)  <  tpiv)  then  jo  to 

if  (jh(i)  -  0)  then  go  to  y 4 4 ; 

if  (kg)  then  go  to  g S 0 ; 

I^O; 

1 4  5 : 

if  (y  (l)  <  dri)  then  go  to  g  S  0 ; 

go  to  g  4  7  ; 

j  4  4  : 

if  (kg)  then  go  to  i4S; 
kg : =  true; 

i  4  7 : 

dd  :  =  y  (i)  ; 
i  r  :  -  i  ; 

1  ‘>0 : 

en  l  g  0  2  ; 

i  f  (ir  *  0)  t hen  <jo  to  l99; 
dd:  =  10. t  2  0 ; 

tor  i:=  1  stef  1  unt  i  1  rr  lo 
1 0  3 :  begin 

if  (y  (i)  <tpiv  V  x(i)s0.0  V  y  (i)  x  ( i)  ) 

t  hen  tjo  to  i  1  U ; 
a.i :  =  x  ( i)  //  (i)  ; 
i  r :  =  i  ; 

]  1  0  :  en'1  j  0  3  ; 

i  f  (  -meg)  then  go  to  j  9  9 ; 
hh:  =  -t.  piv; 

tor  i:=  1  stej  1  until  rr  io 
l04:  beg i n 

i_f  (x  (i)  >0 .  V  y  (i)  •  t  h  V  y(i)*dd>x(i))  t hen  go  to  i  10; 
bn:  =  y  ( i)  ; 
i  r  :  =  i  ; 

] 3  0  :  en  1  l 0 4 ; 

)49:  _i_f  (ir  =  0)  then  to  b07; 

comment  pivot  on  fir,  if)  (step  6)  ; 

i  f  ( ih  (ir)  >0)  then  kt  ( jh  (i  r)  )  :=  0; 
i  0  0 :  Fi vot  (  i  r ,  i  t , e ,  j h  #  k b , pe , x , y ,ir , cost , nun  ov #  ver)  ; 
i_f  (nurrnv  <  rr)  t  hen  go  to  100; 
i  n vc  : -  i  nvc  ♦ 1 ; 

if  (invc  =  nv<>r)  t  h  e  r  i  go  to  rr  2  0 ; 
go  to  100; 


■j’** 
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comment  end  ot  algorithm,  set  exit  values; 
b07:  i_f  (  -»feas  V  rcost<- 1 000 . )  then  cjo  to  b03; 

comment  infinite  solution; 

E  :  =  2 ; 

go  to  b50 ; 

comment  problem  is  cycling  perhaps; 

a  60 ;  k : =  4 ; 

go  to  bSO; 

comment  feasible  or  infeasible  solution; 

b0  3 :  £ ; =  0; 

PSO:  H  (  -tfeas)  t  hen  k:=  k  ♦  1; 

for  j :  =  1  step  1  until  n  do  z  (j)  ;  =  0.0; 

f<  r  i :  =  1  step  1  until  m  TTo  i_f  { jh  (i)  >0) 

then  z(jh(i))  :=  x  (i)  ; 
ko  ( 1)  ;  -=  k ; 
ko  ( 2)  :  =  iter; 
ko  (3)  ;  =  i nvc; 
ko (4)  : =  numvr ; 
ko  (S)  :  =  numpv; 
ko  (6)  :  =  jt ; 
end  SIMPLE; 


procedut e  Get 
begin 

rea 1  ai 
f 00;  lor  i: = 
LL:  =  0; 


Column  (jt ,a, c,e, pe,y , m, ccst , texp, 

jt;  integer  i,j,LL; 

1  step  1  until  m  do.  y  f i )  :=  0.0; 


cost; =  c  (jt)  ; 

for  i :  =  1  step  1  until  ir  do 
f  1 1 :  begin 

aijt;=  a (i, jt)  ; 

i  f  (aijt  =  0.)  then  ^o  to  f02; 
cost:=  cost  ♦  ai  jt  ♦  peTT)  ; 
for  i ;  =  1  step  1  unt  i  1  rr  do 
begin 

LL:  =  IL  ♦  1; 

y  (i)  :  =  y  (j)  ♦  ai  jt  •  e  (LL) 

end ; 

go  to  f  0 5 ; 

f 02 :  LL; =  LL  ♦  m; 

tOS;  end  f  11 ; 

ymax;=  0.0; 

for  i ;  -  1  step  1  unt  il  w  do  (abs  (y  (i)  ) 
then  ymax;=  afcs  (y  (i)  )  ; 
tpi v  ;=  ymax  *  texp; 
eml  Get  Col»imn; 


tp 


iv,  ymax) 


yma  x) 
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rroce  jure  Pivot  (ir,  jt e,  jh ,  kb#  pe,  x,  y,  rr,cost ,  numf v,  ver)  ; 
begin 

red  1  xy,xold,yi;  i  nteger  i,j,L,LL; 
lOO:  nurnpv:=  numpv  ♦  1; 

jh  (ir)  :=  jt; 
kh  ( jt)  :=  ir; 
yi:=  -y  (ir)  ; 
y  (ir)  :  =  -1.0; 

LL:  =  0; 

for  j:=  1  step  1  until  rr  do 
i 0 \ :  Begin 

L:=  LL  ♦  ir; 

if  (e  (L)  /  0.0)  then  go  to  i05; 

•  —  LL  ♦  itf; 
go  to  i()4; 

i  05  :  xy :  =  e  (L)  /yi  ; 

pe(j)  :=  pe(j)  ♦  cost  *  xy; 
e(L)  :=  0.0; 

for  i :  =  1  step  1  unt  i  1  it  do 

beaia 

LL  :=  LL  ♦  1; 

e  (LL)  ;=  e  (LL)  ♦  xy  ♦  y  (i) 

end; 

i04:  end  i 0 1 ; 

xy  ;  =  x  (ir)  /  y i; 
lilt  i :  =  1  step  1  until  m  Jq 
begin 

xold:  =  x  (i)  ; 

x  (i)  :=  xold  ♦  xy  ♦  y  (i)  ; 

if  (  -iver  A  x  (i)  <0.  A  xoldiO.)  then  x  (i)  :=0. 
end; 

y  (ir)  ;  =  -yi  ; 
x  (ir)  :  =  -xy; 
end  Pivot 
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